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Abstract 

The Markov chain Monte Carlo method as a statistical mechanics technique for the study of 
macroscopic systems has furnished the scientific community with great knowledge and advances 
in the theory of phase transitions. While a number of Monte Carlo models have been proposed 
for the study of surface growth, these models have not nearly being studied as exhaustively as 
in the models of magnetic systems, a paradigm of which is the classical model of Ernest Ising. 
In particular, studies of phase transitions in surface/interface science are almost non-existent. 
This article has been written to motivate research in this area of statistical mechanics from the 
perspective of surface science. In this article we survey the rudiments of the method along with 
some models of disordered systems such as magnetic systems, material fracture, nano-pattern 
formation under ion bombardment, and molecular chirality. 
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I. INTRODUCTION 



The relevance and progress of theories of natural phenomena is strongly rooted in exper- 
iments that probe these phenomena under viable and accessible conditions. However, with 
increasing technological advances the experimental apparatus evolve with increasing sophis- 
tication often accompanied by huge operation or acquisition costs that create barriers which 
scientific enquiries must overcome before experiments designed to investigate them can be 
approved. Fortunately, the increasing technological innovations have created an avenue for 
experimental probes through the enablement of numerical simulations that mimic real-life 
experiments in the virtual world embedded in the central processing unit of a computer. 
Such simulations are based on the theory of the observed phenomena and, after verifying 
their results through comparison with experimental results, can be used to study areas not 
yet experimentally studied due to the barriers mentioned above. The results of the sim- 
ulations then determine whether the barriers to the corresponding experiment are real or 
imaginary. 



The Monte Carlo method is a probabilistic simulation met 



lod used in applying the theo- 



ries of statistical physics to the study of macroscopic systems ihsl]. Such systems are termed 
"disordered" due to the large (A^ — ?■ oo) degrees of freedom and the stochastic nature of their 
macroscopic events. Dynamics of surface nanostructure formation and evolution are com- 
monly studied by means of stochastic partial differential equations or discrete solid-on-solid 
simulation models. Each of the sparse (discrete) simulation models are difficult to directly 
map to the continuum space and the statistical mechanics foundation of the solid-on-solid 
models are often left out in their formulation, most of which focus mainly on the calculation 
of the probabilities of the relevant events and the consequent simulation approaches with- 
out much recourse to the statistical mechanics background. This makes it difficult to effect 
the above mapping and leaves a number of unanswered questions. The incompleteness in a 
number of reports, under the general assumption that they are either known or irrelevant, 
occurs in other well. 

However, the statistical mechanics theory is more completely explored and developed 
in the application to magnetic systems and phase transitions Some of these 

advances result of what obtains from the application of the theory to seemingly 

unrelated areas (e.g self-organised criticality in granular mounds, biophysics, etc.), which 



have also benefited from tlie reverse feedback of applied results of the fostered advances in 
statistical theory. 

In this article we review the Markov chain Monte Carlo method in general, as well as 
some models of disordered systems to which it has already been applied. The models consid- 
ered are the Ising model of magnetic systems, a sputter-erosion model for studying pattern 
formation on material surfaces at nanometre length-scales, a model of material failure, and 
a model of molecular matching and packing. These models allow for the application of 
standard techniques of analysis in statistical physics. With recent advances in algorithm de- 
signs and combinatorial optimization in theoretical computer science 0, Q], the simulation 
algorithm for these models are being continually developed with increasingly faster and effi- 
cient simulations [4] which can help develop other fields and vice- versa in the ever-thinning 
interface between the fields of scientific endeavors. 

The rest of this article is presented as follows: In section |TT] we review the basic but 
fundamental features of a computer simulation of disordered systems based on the theory 
of statistical physics. We start from standard textbook material and gradually progress to 
research methods at the level of advanced texts and reference materials; we focus on the 
physics all through and exempt the necessary but avoidable mathematical rigor for ease of 
assimilation and appreciation of the physics. In the context of this established theoretical 
simulation basis, we then discuss a few of the existing models, referred to above, in sections 
nil Al nil Bl nil Cl and IIIIDI to facilitate the creation of new ones in yet unexplored areas. 
And, in particular, to motivate research towards the solution of open problems of statistical 



mechanics of disordered systems 



in surface science (e.g. see Ref. 



rom the perspective of surface science and new frontiers 



II. REVIEW OF THE MARKOV CHAIN MONTE CARLO METHOD 

The macroscopic or thermodynamic properties of a system can be calculated as a function 
of its microscopic properties through the relation (see e.g. Ref. 

S = kB\ndn, (1) 

where the entropy S" is a macroscopic state variable, dQ is the number of microscopic states 
with energy in the interval dE and /c^ is Boltzmann's constant. Note that the statistical 



definition of entropy given in Eq. ([T]) contains an implicit additive constant, and is for 
temperature defined in absolute units. For macroscopic changes of state of the system 
under reversible (infinitesimal) processes the first law of thermodynamics states 



dE = Y,y,dX„ (2) 

i 

where E is the internal energy of the system, i/i are the generalized forces (e.g. intensive 
variables like pressure, surface tension, etc.), and Xi are the corresponding generalized dis- 
placements (e.g. extensive variables like volume, area, etc.) js, 9|. For irreversible processes 
the state variables E, yi, and Xi are undefined. For processes in which the relevant state 
variables are, for instance, S, the volume V, and the number of moles n, then E = E{S, V, n) 
is the function of state and according to Eq. ([2]) the infinitesimal change in the entropy is 
given by 

dS = ^dE + ^dV - ^dn, (3) 
to which an application of Euler's relation gives 

5 = 1 + ^^^-^. (4) 

Depending on the nature of the interaction of the system with its environment there 
may be one or more of energy exchange, matter exchange, or exchange of spatial extent, 
which determines the statistical ensemble, i.e. assembly of microstates subject to one or 
more of these constraints, within which the system evolves and is simulated. In any such 
scenario the equilibrium state of the system is any of a large number of equally probable 
microstates, whereas the probability of a microstate varies between equilibrium states. In 
fact the number of microstates is infinite if the phase space points form a continuum. 

For reversible processes carried out on a canonical ensemble where V and the number of 
particles A^(A^ oc n) are constant, dV = dN = and the number of microstates Qi with 
energy £i at temperature T is given by 

= exp ■ (5) 



The relationship between E and Si is that between the average of a data set and the 
elements of the data set; as we shall see below the system must be allowed to equilibrate at 



a particular temperature before measurements are made, nevertheless, the energy fluctuates 
about an average value after the equilibration time such that E = J2j ^j] where u is the 
number of microstates with (possibly different) energies £j through which the system evolved 
in the measurement process. According to the fundamental postulate of statistical mechanics 
all the microstates corresponding to a particular macroscopic state of an isolated system 
are equally probable (microcanonical distribution law), and at small times the canonical 
ensemble behaves as if microcanonical, thus the probability Pi of a microstate with energy 
Si, Pi oc 1/Qi, i.e. 

P. oc exp (-^) , (6) 

or 

where the constant Z that ensures the normalization of the probabilities Pi, so that 
X^j -Pj = 1, is the so-called partition function given by 

^-i:exp(-^). (8) 

Note that Q = constant for the micro-canonical ensemble, which implies Pi = constant 
for this ensemble, hence, such states can be sampled uniformly as in Monte Carlo integration 
and not by importance sampling as is the case in the canonical ensemble. 

The requirement that the macroscopic states be reversible imposes the condition of er- 
godicity on the microscopic phase space of this system so that every microstate must be 
accessible along a flnite path or trajectory that connects the microstates in phase space. 
This is achieved by constraining the simulated microstates to a Markov chain such that a 
new state is attained from a previous one in the chain. Thus in simulations of disordered 
systems in a given statistical ensemble the goal is to visit the phase points or ensemble 
microstates m with a probability Pi proportional to some given distribution p(m), where 
the phase points are locally correlated in a sequence of points mi, m2, m3, . . . each derived 
from the subsequent one and all constituting a Markov chain. 

In the canonical ensemble p(m) is given by the Gibbs/canonical distribution function 
A; exp (^—-^j^^ and Pi is given by Eq. [71 p(m) is not necessarily normalized to have unity 

5 



integral or sum over the sampled region but is proportional to a probability. Now Pi is a 
conditional probability, i.e. the probability of accessing the microstate rrij given that the 
system is in the microstate mj_i 



Pi = P(mi|mi_i) oc p(mj) (9) 

So Pi are transition probabilities of transition from microstate mj_i to rrij . For us to 
obtain the partition function we would have to first sample all microstates and evaluate 
the sum of the distribution functions p(nij) which is, however, unnecessary since (assuming 
1 or by simply using the normalized distribution) we must have 

P(m^|m,_i) ^ pinii)/Z 
P(mi_i|mi) p(mi_i)/Z 

that is, 

p(mi_i)P(mi|mi_i) = p(mj)P(mi_i|mi) (11) 

The condition ( ITTl) is known as detailed balance. It allows us to calculate transition prob- 
abilities in terms of relative transition probabilities between successive microstates without 
actually calculating the partition function. For instance, in the canonical ensemble Eqs. 
and ffTOj) imply 

Pm,.,-.m, = exp (-^) . (12) 



where AS = £i—£i-i is the energy difference in a transition from the microstate mj_i with 
energy £i^i to the next microstate nij with energy £i in the Markov chain, and Pmi_i^mi = 
p(mi|mi_i) ^j^g relative conditional probability that the system will evolve to microstate 
nij from the microstate mj_i. 

According to Eq. f lT2|) the system will certainly change into the microstate irij from 
mj_i if A£i < A£i_i since in this case Pm^.i-i-mi > !> whereas it will likely remain in the 
microstate mj_i if A£i > A£i_i since Pm,_i^m, < 1 in this case and Pmi_i^mi — > as 
A£ — i- oo. This ensures that the vast phase space is "importance" sampled in such a way 
that the system evolves in the few most probable microstates in the phase space in the course 
of the simulation and reduces the rigor and inefficiency of having to go over or uniformly 



sample ("simple" sampling) the entire phase space in the calculation of average quantities 
such as, e.g. 



. y-^-p (13) 

In the calculation of Eq. fll3p most of the Gibbs weights exp {—(5£i) are vanishingly small, 
corresponding to the very large number of microstates in which the system is unlikely to be 
found, and therefore contribute very little to the average value. Note that, depending on 
the application, £ and T may have entirely different physical meaning. 



A. Simulation Algorithm 

Based on this theory, a Markov chain MC simulation of the evolution of the microstates of 
a system in statistical mechanics can be performed by using the following algorithm, known 
as the Metropolis algorithm: Generate a trial configuration mj_i 



m,_i = {(Ti,(T2, . . . ,(TAr} (14) 

for the initial microstate of the A^-particle system and obtain the next microstate rrij in 
the Markov chain by tweaking the coordinate ak of a single ( A;th) particle. Calculate the 
acceptance probability a;(mj_i, mj) of the new microstate rrij using the formula 



a(mi_i, rwi) = min (l, Pm,_i^m,) (15) 

This means that the transition mj_i — >■ irij is accepted if Si < Si-i [a(mj_i,mj) = 1] 
otherwise it is accepted with a probability [a(mj_i, rrij) = Pm.^i^mi < 1] which, at constant 
T, decreases with increasing AS. 

This algorithm therefore ensures that the system always tend towards microstates of 
lower energy than the previous one but unlike the so-called "greedy algorithms" 10|, that 
never accept a higher energy option, it also allow the system to sometimes (though with 
very low probability) accept a move to a higher energy microstate; wherein lies its power 
and vast applications to diverse optimization (minimization/maximization) problems. For 
instance, if higher energy moves are always rejected then there will be no possibility for the 
system to move out of any of the local minima in Fig. [1] and the system would therefore not 



equilibrate to the global minimum. However, as the temperature decreases the probability 
of the system moving out of a local minimum, if trapped there at low T, also decrease and at 
such low T it may be difficult to attain a microstate of global energy minimum; depending 
on the energy landscape of the system. 

It is best to run independent simulations with different initial conditions, e.g. one with 
a cold start (T = 0) and another with a hot start (T = oo), so as to avoid inaccurate 
results due to local minima pitfalls. For the cold start a starting configuration or microstate 
that is known to have the lowest energy is used, while for the hot start a configuration that 
corresponds to the highest energy is used. The result of the independent simulations must be 
the same irrespective of the starting configuration, otherwise one or both must have ended 
in local minimum and further checks are necessary. 




FIG. 1. Hypothetical energy landscape of two systems. In (b) the global minimum is within a 
very narrow region of the phase space of the system which means that reaching it may be almost 
impossible if trapped elsewhere at low . 
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Since the system equilibrates to a "global" energy minimum in which it then spends the 
rest of its time, the probability Pm,^m,_i is time dependent , i.e. Pmi-i>mi_i = Pmiit) and 
its time dependence is governed by the master equation 



- [-Pm, (i)-Rm,^m, - -Pm, (^)-Rm,^mJ (16) 

where Pmi{t) is the probability of the system being in the microstate rrij at time t, and 

dP (t) 

Raii^mj is the rate of transition from microstate mj to uij. In equilibrium — ^ — = and 

Obviously the efficiency of this algorithm, at the so-called "importance sampling" of 
the phase space region in which the system spends most of its time, will depend on the 
energy landscape of the system; the algorithm will be more efficient and faster for a system 
with the energy landscape of Fig. [T]^a) than that of Fig. [U^b). As a result, a number 
of cluster algorithms have been proposed to increase the simulation speed and efficiency at 
"importance sampling" by avoiding the slow changes of state that occur with time in certain 
stages of the simulation (e.g. critical slowing down, frustration). Such algorithms consider 
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M- 



microstate transitions involving clusters of particles instead of single particles 

Since the system may end up in local minima it is better to let the simulation run for a 
while before taking measurements so as to allow the system to equilibrate; the time required 
for the system to equilibrate is known as the equilibration time, which depends on the energy 
landscape of the system. The configurations, or microstates, are usually stored as arrays 
whose elements are the single-particle states corresponding to each particle of the system. 
The energy of a microstate depends on the interactions between the particles and is the sum 
of the energies of these interactions. 

The particles, whether distinguishable (e.g. Maxwell-Boltzmann particles) or indistin- 
guishable (e.g. fermions or bosons), must have identical surroundings and interactions 
in such a way that the boundary particles are not distinct from the inner ones. This is 
achieved by imposing periodic boundary conditions on the system, thus transforming the 
D— dimensional system into a {D + 1)— dimensional torus, so that the particles at one end 
of the system sees the particles at the other end as nearest neighbors. The storage require- 
ments for specific problems or the geometrical nature of certain systems may require other 
boundary conditions like screw-periodic, anti-periodic, mean-field, or free-edge etc, bound- 



ary conditions. Other lattice structures may also be more appropriate for specific problems 



a 



B. Uniform Deviates 



The recipe for simulation is incomplete without (pseudo) random numbers 10|. We can 
not choose at in a deterministic version because that would be unrealistic, it would either 
imply a prior knowledge of the phase trajectory followed by nature or that we are the 
ones dictating what course nature would take which is contradictory. Hence, is chosen 
randomly by using floating-point uniform deviates such that each at has the same likelihood 
of being picked as any other. Our knowledge of the probability of a trial microstate is also 
not enough because we need to know when to accept or reject it according to its probability. 
Again, a floating-point uniform deviate is used. 

For instance, if the probability of an event is 0.1 then the event will occur once out of 
every 10 attempts, as the number of attempts tends to infinity. To simulate this either 
of two approaches may be used, in both approaches the random numbers are assumed to 
be uniformly spread over an interval (0 to 1 for floating-point random numbers) such that 
each random number has an equal probability of being chosen as any other one. In the 
first approach the interval is divided into ten shorter intervals of equal width (e.g. — 0.1, 
0.1 — 0.2, 0.2 — 0.3, etc; either excluding lower bounds or excluding upper bounds), the 
event is associated with one of the intervals (corresponds to one face of a ten-faced die), 
random numbers are generated one at a time (analogous to throwing a ten-faced die), and 
the event is taken to occur when the random number falls within the associated interval. If 
the probability is e.g. 0.7 then the occurrence will be associated with seven intervals out of 
ten or the non-occurrence may be associated with three intervals. In the second approach 
the random numbers are associated with either an occurrence or a non-occurrence of the 
event. If the random number is less than or equal to the probability of the event (e.g. 0.1) 
the event is said to have occurred else the event has not occurred. The later is the one 
adopted in simulations, and is the same as the former since their acceptance rates are the 
same. 
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C. Finite Size Effects 



The number of particles (system sizes) that may be considered in a simulation is much 
less than the number of particles in a typical macroscopic system, i.e. the simulated system 
size << oo whereas in the thermodjTiamic limit we should have N — t- oo. Typically, 
real ~ 10^^ particles whereas simulation ^ 10^. For instance, a system with only two 
possible single-particle states will have 2^ ^ 2^°^^ microstates or configurations while the 
simulation phase space comprises of ~ 2^*^° microstates. 

However, use of small N still gives very accurate results when compared with real systems 
for which N ^ oo provided the simulation is for temperatures far from the critical temper- 
ature at which a phase transition occurs. For temperatures close to and above the critical 
temperature Tc the effects of using small A^, known as finite size effects, become very strong 
and the simulation results become inaccurate; the inaccuracy increases with decreasing A^. 
Hence, for studies involving phase transitions, or for simulations at temperatures above Tc, 
it is important to do a finite size scaling of the results to obtain the actual results of a 
real thermodynamic system. The form of the finite size scaling will depend on the set of 



thermodynamic quantities being measured and their scaling relations [l-3|. It also depends 
on the order parameter, which is a quantity that has a value of zero on one side of the phase 
transition and a nonzero value on the other side of the phase transition; i.e. which is like a 
step function whose step occurs at the phase transition. 

The origin of the finite size effects is in the divergence of the correlation length ^ (the 
length scale of the fluctuations in the order parameter) at the critical temperature as a 
result of which it is impossible to get accurate results of the thermodynamic properties of 
the system at from any simulation performed on a finite lattice. Since the lattice size is 
less than the value of C, at Tc, all the lattice points are correlated and one is unable to obtain 
the scenario of linearly independent single-particle states. 

III. MODELS OF SOME DISORDERED SYSTEMS 
A. Magnetism and the Ising Model 



Magnetism is one of the most studied phenomena in materials. It has been a study 
of active and intense research interest for centuries jisl . flG^ . Any material is composed 
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of atoms bound together by cohesive van-der-Waals, covalent, electrovalent, or metaUic 
(electron gas) interactions which reduce to electrostatic interactions between dipoles made 
up of the negatively charged electrons on one end and the positively charged ion cores on 
the other. Each electron has an intrinsic, purely quantum mechanical, angular momentum 
called the spin S of the electron. Associated with each electron spin is a magnetic moment 
of magnitude (i.e. noting that the gyro-magnetic ratio ^ in this case is twice the 
ratio of the magnetic moment to the angular momentum) which is solely responsible for 
the magnetism of the material; e is the electronic charge, m is the mass of the electron, 
and c the speed of light in vacuum. Thus the magnetic property is a quantum mechanical 
property and docs not have a classical analog, which means that the Hamiltonian "H of the 
system does not include position and linear momentum operators. The exchange energy due 

2 

to spin-spin interactions between spin i and spin j can be specified as where rij is 

the distance between the spins. For a magnetic material in a magnetic field, therefore, the 
only changes in the energy of the system are those due to the interaction of the spins with 
each other and with the external magnetic field H. Hence, the Hamiltonian of the system 
is given by 

^ = ^ E E -Si^^ - — Y.HSi (17) 

This is equivalent to the model of Ernest Ising for which, for simulations on lattices of unit 
spacing with nearest neighbor interactions, the eigenvalues of the Hamiltonian is expressed 
as 

E^-JY.S,S^-HY,S,, (18) 

{hi) ' 

where J is the unit of the spin-spin exchange interaction and is often chosen to be -|-1, 
in which case the lowest energy state will tend to have neighbouring spins aligned and the 
low temperature ground state will be a ferromagnet (all spins aligned), or —1 in which case 
the lowest energy state will tend to have neighboring spins opposed in orientation and the 
low temperature ground state will be an anti-ferromagnet (alternating spins), (i, j) denotes 
a nearest neighbor pair, and H is the unit of the spin-field coupling. 

The spin operators (actually S — Sz for each spin, assuming the field is oriented along 
the z-axis so that H — Hz) have eigenvalues where h is the Dirac action constant. In 
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simulations the eigenvalues of the spin operator 5* are chosen to be ±1. Thus the number of 
single-particle microstates of the system is reduced to only two; for up and down orientations 
of the intrinsic spin. When the exchange interaction energy J varies from bond to bond then 
we have a spin-glass system which is widely studied with a variety of interesting ramifications 
17|. 



B. Surface Nanostructure Formation and the HKGK Model 
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FIG. 2. Surface profiles for different sputtering times (Monte Carlo steps). From top to bottom, 
left to right: t = 0.5, 1.5, 4, 9, 14, and 20 ions/atom. The vertical scale indicates the surface height 
range. The ion beam direction is indicated by the bar on the profiles, showing that the ripple 
orientation is perpendicular to ion beam direction as expected since 9 = 50° here (After Ref. [18]). 
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A more recent observed phenomenon, in surface science, is the formation of self-organized 
nano-patterns on material surfaces when bombarded by a beam of energetic particles, with 
energies in the range of keV 0, Bfl- 

The surface morphology arises as a result of the 
interplay between the roughening processes of stochastic removal of surface particles by 
ion-sputtering, and the smoothening processes of diffusion of surface particles or radiation- 
induced viscous flow. The roughening process creates instability of the surface against 
further perturbations by the eroding ions such that the sputter yield becomes curvature 
dependent with the result that surface depressions are eroded in preference to surface pro- 
trusions. 

In the continuum theory the surface evolution is modelled with deterministic and stochas- 



tic partial differential equations 



22 



-24 



26| . According to the continuum theory, it is the 



curvature dependence that leads to the formation of ripples of wavelength A = 27ry 2i^/|z/|, 
where K is the surface diffusivity, and z/ is a surface tension coefficient. The wavelength 
of the ripples is on the nanometer scale and they are oriented along the direction with the 
largest \v\ 22]. The ripple orientation is parallel to the projection of the ion beam direction 
onto the surface plane for large incidence angles (close to grazing incidence), and perpendic- 
ular to the projection for small incidence angles; except for metallic surface with anisotropic 
diffusion in which case the ripple orientation is perpendicular to a crystallographic direction 
(i.e. the one favored for diffusion) at small incidence angle. 

Hartmann, Kree, Geyer, and Koelbel (HKGK) proposed a discrete solid-on-solid Monte 
Carlo model 27|] for the simulation of surface modification and evolution by ion-bombardment 
which, in contrast to the Ising model above, simulates the surface morphology indirectly 
by resorting to the effect of collision cascades on the surface evolution without the direct 
need for a surface Hamiltonian. These collision cascades are set off in the near sub-surface 
layer by the impinging ion and are responsible for the redistribution of the energy of the 
impinging ion to the surface particles. The energy -E(x) received by a surface particle at 
position X = (xi, X2, x^) due to the arrival of the ion somewhere in the vicinity of this surface 
particle is assumed to be of the Gaussian form [28|: 



where a and p are the widths of the collision cascade ellipsoid parallel and perpendicular 
to the ion beam direction, respectively. A sputtering process is simulated by eroding surface 



X2 I ^2 
1 "I" 
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^articles in the vicinity of the impact position of the ion with probabihties proportional to 
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29 



30| . and a diffusion or smoothing process is simulated by using any of the available 



diffusion models [e.g. Refs. 3l|, |32|], depending on the experimental hopping rates for the 



material concerned. Profiles obtained from such a simulation are shown in Fig. |5J 



C. Molecular Chirality and the YNL Matching Algorithm 



The problem o 
its mirror image 



measuring the extent to which an object fails to coincide with or match 
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- |41| when optimally superimposed on the mirror image was recently 



tackled in a direct way through simulations of the matching process using the Markov chain 



MC method 42| (YNL algorithm). The measure of chirality or degree of matching of the 
objects (an object and its image or two different objects as in molecular similarity or facial 
resemblance) that have been optimally superimposed by the YNL algorithm is based on the 



34| 



concept of the Hausdorff distance between sets as proposed by Buda and Mislow 

This direct approach entails the simulation of the matching process which includes a large 
number of rotations and translations such that the Hausdorff distance between the objects 
being matched is as small as possible, dn = if an object is achiral or if two objects match, 
and dn > if an object is chiral or two objects do not match. The challenge is that one may 
get very close to the global minimum of dn only to be pushed into a completely different 
region of the dn landscape with the slightest of translations in the matching process. But 



the YNL Monte Carlo matching algorithm 



42| always arrive at the optimal match. 



The YNL algorithm searches for the optimal superimposition of an object A and its 
mirror image A' = —A, or second object, by performing a series of random rotations and 
translations of A keeping A' stationary. The simulation is started with a random orientation 
and translation of A which represents the initial state or configuration. The matching process 
is iterated down over the 'temperature', thus annealing the process and ensuring that as the 
system cools, the rotation and translation reduces in magnitude so that at high 'temperature' 
the system moves freely from one region of the d^ landscape to another in the search for the 
region that contains the global minimum, and much less freely when this region is supposed 
to have been found (at low temperature). During each iteration stage over the temperature, 
it is kept fixed while Nmc Monte Carlo steps of the matching process are performed. 

In each MC step the state of the system, which is the orientation and position of A 
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relative to A', is changed by rotating and translating A and the move is accepted according 
to the standard MC transition acceptance rules, tailored towards acquiring minimum dn 
or dff^.^ provides a measure of the chirality xh of the object through the chirality 



measure of Buda et al 



35j: 



dn ■ (A, A') , , 

■,„(A) = J!^j^ (20) 

where du^:^^ is the minimum Hausdorff distance between A and A over all positions 

and orientations of A and A', (i(y4) is the diameter of A. 



D. Material Fracture and the Fibre Bundle Model 



Fracture in heterogeneous materia 
of research interest for a long time 



s is a complex physical problem which has been a focus 



43l-l45l|. Fibre bundle models (FBMs) form a fundamental 
class of approaches to the fracture problem through their capture of the essential features 
of material breakdown. They are models of materials as composites of fibres such that the 
breakdown of the material is as a result of the propagation of fracture among the failing 
fibres. 

One such model which provides a deep understanding of the intrinsic nature of the fracture 
process is the dynamic fibre bundle model According to the dynamic FBM the material 
breaks down as a result of the fatigue of its fibres over time. That is, when the material 
is stressed or loaded the fibres share the load and get progressively stressed with time such 
that even if the external load or source of stress is removed the acquired fatigue in the 
fibres remain and are incremented with the introduction of a new source of stress until the 
individual fibres reach their threshold and f ail, leading to the overall failure of the material 



47| which attributes the failure of the individual 



with time. This differs from the static FBM 
fibres to their quasi-static loading due to the gradual increase in the external load. 
The dynamic FBM includes a normalised combined load sharing rule j46l] 



S^^^^n) = ^^A—, (21) 



Tij is the distance between an active fibre i and a failed fibre j, F denotes the set of active 
fibres, and 7 is a variable parameter that controls the effective range of interaction among 
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the fibres. This combines the local load sharing rule (7 = 00) in which the load borne by a 
failed fibre is shared equally among the four nearest neighbour fibres, the global load sharing 
rule (7 = 0) in which the load of a failed fibre is shared among all existing fibres, and the 
variable range load sharing rule (0 < 7 < 00) in which the load of a failed fibre is shared 
beyond the nearest neighbours but not globally, with the range depending on 7. 

The probabihty Pj{t) of the breaking of a fibre j in one sweep of the lattice in the time 
interval Sk is given by 

P,{t) = (22) 

p is the WeibuU index (2 < p < 50) which gives the degree of heterogeneity of the system; 
as p increases the material becomes more homogeneous. crj{t) is the stress on the fibre j at 
time t. 

IV. CONCLUSION 

This article reviews the modern theory and implementation of computer simulations of 
disordered systems, based on the statistical mechanics approximation method of Markov 
chain Monte Carlo techniques. We started with the fundamentals without recourse to any 
system or problem, for a generalized treatment, and then discussed a number of recent 
applications including Ising model or spin-glass model of magnetic systems, a solid-on-solid 
model of nanostructure formation on surfaces driven by ion-bombardment, statistical MC 
models of material fracture, and a MC model of chirality which is applicable as a matching 
algorithm. 
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